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1 Introduction 

Modern methods for rare events allow the search of metastable states in complex 
systems, such as systems in condensed phase. Examples of these methods are the 
Temperature Accelerated Molecular Dynamics (TAMD) (TJ, with the associated 
Temperature Accelerate Monte Carlo (TAMC) |2] and the Methadynamics EE). 
These methods are based on the introduction of an extended system consisting of 
the physical degrees of freedom (positions and momenta of the atoms) and a set of 
extra dynamical variables z = {z/}/=i, m representing possible realizations of a set 
of suitable collective variables {ft(x)}/ = i m characterizing the states of the system. 
From now on, we will refer to the space of the z variables as the z-sp ace. The z 
variables are coupled to the atoms by a suitable potential (see Sec. |2 1| ) that forces 
the latter to stay in a configuration compatible with the restraint |0/(x) ~ zt}i. 
From these simulations we can compute the probability density function P(z) to 
be at a given point in the z-space. P(z) is connected to the Landau free energy by 
the relation F(z) = — &#riogP(z), where ks is the Boltzmann constant. In these 
methods, the sampling of the configuration space is accelerated by biasing the dy- 
namics of the z variables: the latter are forced to visit less probable regions of the 
z-space and, as a results, the atoms will also visit the corresponding regions in 
the configuration space. In TAMD and TAMC this biasing consists in setting the 
z- temperature to a value much higher than the physical temperature. This makes 
it possible for the system to sample regions at high free energy, F(z). In Metady- 
namics, the biasing amounts to apply a history dependent biasing potential on the 
z variables that forces the system to visit region of the z- space not visited before. 
In both cases, the z variables, and as a result the atoms, can more easily overcome 
free energy barriers present in the system. 

Despite this acceleration, the reconstruction of the free energy surface F(z) 
by "binning and histogramming" is very expensive. However, often we are in- 
terested only in the identification of metastable states and the characterization of 
their relative stability. This task can be accomplished by following a more effi- 
cient, two steps approach. In the first step the z-space is scanned by relatively 
short TAMD/TAMC or Metadynamics runs. The corresponding z-space trajectory 
is analyzed to find "clusters" 0[6]|, i.e. those regions in the z-space where the tra- 
jectory variables tends to spend a significant amount of time. These runs must be 
long enough to allow the identification of metastable states on F(z) but can be sig- 
nificantly shorter than the length required to accurately reconstruct the free energy 
surface. In the second step, the difference in free energy of two of the metastable 
states identified in the previous step is calculated by numerically integrating the 
mean force VF(z) along a path connecting them (thermodynamic integration). 
By repeating this operation for a proper set of pairs of states it is possible to re- 
construct the relative free energy of all the metastable states. The VF(z) can be 
estimated by computing the ensemble average conditional to {ft-(x) = z/}/ of a 
suitable observable (see Sec. [2]). This task can be performed by restrained (con- 
strained (71) MD or MC (llEl. This approach is based on the assumption that there 
are no free energy barriers in the domain orthogonal to the z-space, as otherwise 
the sampling of the conditional probability density function p(x|{0j(x) = z/}/) 
would be quite inefficient. We explain the problem of sampling p(x|{0j(x) = z/}/) 
by the following example. Consider a free energy surface of the type shown in 
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Fig. 1 (Left) Free energy landscape of a system with two slow degrees of freedom represented in 
colormap: blue — » low free energy, red — » high free energy. (Right) Free energy surface projected 
along Z2 at 6\ (x) = z\ 

Fig. [I] In this system there are two slow degrees of freedom, 0(x)i and 0(x)2, but 
the two metastable states can be characterized by only (x) i : the first corresponds 
to the condition 0(x)i = z[, and the second one corresponds to 0(x)i = z f (. This 
means that, indeed, we can determine their relative stability by integrating VF(zi) 
between z[ and z![. However, we notice that the free energy profile in the direction 
Z2 for a fixed value of z\ (F{zi | 0i (x) —z\)) presents a barrier higher than the ther- 
mal energy. As a result, the sampling of p(x|0i (x) = z\) by restrained/constrained 
MD, or any other standard method, is biased. 

In principle, this problem could be solved by including all the slow degrees of 
freedom in the z- space. However, already the identification of the collective vari- 
ables characterizing the metastable states is very difficult and the identification 
of all the slow degrees of freedoms along the paths connecting these states often 
results just impossible. In this paper we illustrate how it is possible, by combining 
restrained MD with the parallel tempering l8l! MT0l[TTl[T2l . to compute the differ- 
ence in free energy between two metastable states even in presence of unknown 
slow degrees of freedom. We apply this method to the study of the phase diagram 
of Si nanoparticles embedded in amorphous silica (a-Si02). We study the nature 
of the most stable phase as a function of the size and temperature of the nanopar- 
ticle. Anticipating our results, we found that, at variance with bulk Si, in small 
nanoparticles (0 < 2.6 nm) at low T (T < 1500 K) the most stable phase is the 
disordered one. For larger nanoparticles the bulk-like behavior is recovered. 

The paper is organized as follows. In Sec. [2] we introduce the combined re- 
strained MD/Parallel Tempering method (RMD-TP). In Sec. [3] we presents the 
results of the study of the phase diagram of the Si nanoparticles. Finally, in Sec. [5] 
we draw some conclusions. 



2 Theoretical Background 

This section is divided into two subsections. In the first one (Sec. |2.1| ) we revise 
the restrained MD method and demonstrate that, under the hypothesis that there 
are no slow degrees of freedom apart those listed in the set of collective variables 
{0;(x)};, this method allows to accurately and efficiently estimate the mean force 
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VF(z). In the second one (Sec. 2.2) we quickly summarize the parallel temper- 
ing method and explain how, by combining it with restrained MD, it is possible 
to efficiently compute the mean force VF(z) even in presence of unknown slow 
degrees of freedom. 



2.1 Restrained MD 

The restrained MD method considered in this paper is a specialization to the case 
of fixed z of the TAMD method of Maragliano and Vanden-Eijnden d. Let us 
introduce the following equation of motion for the atoms: 

mx = -V x U k (x,z) +thermo(p) (1) 

where m is the physical mass, /3 = (&#r) _1 and thermo(P) indicates that the 
atoms are coupled to a thermostat at /3. The potential U k (x, z) = V (x) +Y4L1 k/2(6i(x) — 
Zi) 2 is the sum of the physical and the restraining potential. Let us compute the fol- 
lowing average: 



m = Hm i ft, »««.»-») = M^M^gfcffl^ 

where 3f k (z) — J dxexp[— j3U k (x : z)]. The second equality in Eq. [2] stems from 
the assumption that, apart for the z, the system is ergodic. f k (z) is the z-th compo- 
nent of the gradient of F k (z) = -jS" 1 log[ J*(z) / JT], where 2T = J dxexp[-/3V(x)] 
is the canonical partition function of the real system. Since is z-independent its 
introduction does not affect our argument but it is necessary for the interpretation 
of Ffc(z) as a free energy. Noting that lim^ exp[-j3 \ ( 0;(x) - Zi ) 2 } / y/{2n/fik) = 
S(0i(x) - Zi), in the limit of k -> 00 F k (z) -/3" 1 log[P(z)] = F(z). In conclu- 
sion, the time average of Eq.|2]along the restrained MD of Eq.[T]is an estimate of 
VF(z), i.e. the gradient of the free energy at the point z. 

Let us now discuss the case in which not all the slow collective variables have 
been identified. For sake of clarity, consider the case in which there are two col- 
lective variables, 0i(x) and O2 (x), but only the first is considered in the restrained 
MD calculation. We further assume that, even if 0\ (x) is not the only slow col- 
lective variable, it is a good order parameter for the process under investigation. 
In this case the free energy profile might look like those reported in the top pan- 
els of Fig. [2] Let us focus first on the one shown in Fig. [2|A and let us compute 
F{z\) by performing the thermodynamics integration on only the VF(zi) along 
z\ as computed by restrained MD. If we start the calculation of the VF{z\) from 
the minimum at low zi, because of the presence of a free energy barrier along zi, 
we shall be sampling only along the valley indicated by the dashed path. We keep 
sampling along this valley up to the point in which the barrier separating the two 
valleys is lower than the thermal energy. At this point the system will collapse 
in the other basin and, by increasing z\ , it will go toward the second minimum. 
Pushing the system forward along z\ it will start climbing the free energy surface 
from the new basin. In short, we would compute the VF(zi) only by sampling 
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Fig. 2 (Top panels) Two cases of free energy landscape of a system with two slow degrees of 
freedom represented in colormap: blue — > low free energy, red — >> high free energy. The dashed 
(dotted) curve represents the part of the space sampled by a pure restrained MD simulation going 
in the direction low-to-high (high-to-low) z\ . (Bottom panels) Corresponding free energy curves 
F(z\) as obtained from thermodynamics integration of the mean force computed by restrain 
MD. 



p(x|0i(x) = z\) around the dashed path sketched in Fig. [2jA. The VF{z\) com- 
puted along the segment 1-2 of the dashed path will be negative and, correspond- 
ingly, the free energy will be decreasing. In the segment 2-3 the VF{z\) will be 
positive and the free energy increasing. This give raise to the minimum at low Z\ . 
In the segment 3-4 and 4-5 VF(z\) will be negative and positive, respectively, and 
the free energy will present a second minimum. The free energy profile associated 
to this path is illustrated in the bottom panel of Fig. [2|A and is represented by a 
dashed line with an arrow indicating the direction along which the thermodynamic 
integration is performed. Let us not perform the same ideal experiment going in 
the other direction: from high to low z\ . In this case the path, denoted by a dot- 
ted line, will follow the other valley. As a result the free energy barrier is shifted 
even if the position of the minima is correctly computed. However, the relative 
free energy of the two metastable states is wrong following both the low-to-high 
and the high-to-low paths as the VF(zi) was obtained from a bad sampling of 
p(x|0i(x) = z\). However, we must consider also the case, shown in Fig.[2jB, in 
which the two valleys are very long and the point in which the system collapse 
in the other basin (segment 3-4) is very close, or beyond, the second minimum 
(see dashed and dotted path in panel B). In the first case the mean force will be 
negative but small along this segment, while in the second case it will be positive. 
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Depending on the exact placement of this segment, the reconstructed free energy 
profile might present one deep and one shallow minimum or just one minimum, 
as in the example of Fig.(2|B. This second example illustrates that also the quali- 
tative picture of the free energy surface can be falsified if not all the slow degrees 
of freedom are taken into account, despite the ones considered could be the good 
order parameters for the process considered. In the next section, we will illustrate 
how this problem can be solved by combining restrained MD with the Parallel 
Tempering. 



2.2 Parallel Tempering and combined Restrained Parallel Tempering MD 

In parallel tempering the simulated system is composed of M replicas of the orig- 
inal system. Each replica is in the canonical ensemble at a temperature 7} (i is the 
index of the replica). Since the replicas do not interact with each other, the proba- 
bility density function that the system is in the configuration xi in the replica 1, in 
the configuration X2 in the replica 2, etc. is: 



w( Xl ,...,x M )- mTu ^ Tu) (2) 

where V (x) is the potential governing the replicas and $ = l/foTJ. £t(T\ , . . . , Tm) = 
YlfLi J dx t exp[— PiV (xi)] is the overall configuration integral. To sample w(x\ , . . . , x M ) 
we can use a MC procedure with two types of moves: (i) the standard MC moves 
within each replica (the specific type of moves depend on the system: atomic, 
molecular, etc.) and (ii) a move in which the positions of the atoms of two replicas 
are swapped. The first type of moves is accepted/rejected according to the usual 
Metropolis criterion: 

a(xi -+ xj) = min ( 1, w{ * h ' ' ' ) = min (l,exp[-ft(V«) - V(x i ))p) 

\ W(X1,...,X/,...,X M J/ 

The second type of moves, which is attempted from time to time during the 
simulation, is accepted according to the probability 



/ w(xi,...,x 7 ,...,x/,...,^ 

a(Xi —> x/,x/ —> Xi) — min 1, — 

V w(xi,...,x/,...,x y ,...,} 



mm(l,exp[-(/3 7 --A)(V(x,-)-V(x 7 -))]) (4) 



where i and j are the indexes of the replicas for which the swap is attempted. Let 
us assume that 7) » 7}. Then, the replica j has an higher probability to overcome 
(free) energy barriers during the standard MC part of the simulation. Therefore, 
it has an higher probability to visit the many metastable states possibly present in 
the system. The swapping move, then, allows the low T replica to sample these 
states too. This explains why, when the system presents energy barriers, parallel 
tempering with M replicas is more efficient than standard MC of length M x Npt, 
where Npt is the number of steps of the parallel tempering simulation. 

Parallel tempering can also be implemented using MD. In this case, the stan- 
dard MC part of the simulation is replaced by a constant temperature MD. As in 
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the MC based parallel tempering, the MD simulations of the replicas are carried 
out in parallel for a given number of steps after which a swapping is attempted. 
This move is accepted/rejected according to the same a(x t — >> Xj,Xj —> x ; ) prob- 
ability reported in Eq. [4] In the MD based parallel tempering, after an accepted 
swapping, the momenta of the atoms are scaled by a factor \]T n ewlT id' Finally, 
both for MD and MC, after a swapping it is necessary to wait for the relax- 
ation of the system at the new temperature before using its trajectory to sample 
w(xi,...,x M ). 

In this work, we want to sample the probability density function of the sys- 
tem conditional to {ft-(x) = z*}*=i, m (see previous section). Consistently with the 
restrained MD method explained in Sec. 2.1 this can be obtained by replacing 
the potential V(x) in Eq.|4|with the potential U k (x,z) introduced in the previous 
section. For an adequate selection of the number of replicas and their tempera- 
tures, the piece- wise RMD-PT trajectories will sample the conditional probability 
density function p(x|{0j(x) = Zi}i) also in presence of additional slow degrees of 
freedom. 



3 Application of the RMD-PT approach to the study of the phase diagram of 
embedded silicon nanoparticles. 

In this section we present the results of our RMD-PT simulations aimed at study- 
ing the phase diagram of a Si nanoparticle embedded in a-Si02. By integrating the 
mean force computed according to the RMD-PT method along a path connecting 
the ordered and disordered states of a Si nanoparticle of a given size and at a 
given temperature, we are able to identify which is the most stable phase in the 
given conditions. The implementation of this approach requires the introduction 
of two collective variables, one controlling the size of the nanoparticle and the 
other controlling its degree of order. These collective variables are introduced in 



the Sec s.|3.1| and |3.2| In Sec. |3.3| we describe the setup of our simulations. Finally, 
in Sec. |4j we present our results, and analyze and discuss the effect of the parallel 
tempering. 



3.1 The M(x) collective variable to monitor the size of the nanoparticle 

We introduce the notion of size of the nanoparticle by a collective variable denoted 
by the symbol £%(x). In a system in which the nanoparticle is made of atoms of 
type 'A' (Si in this case) and the matrix is made, or contains, atoms of type 'B' 
(O in this case) a possible definition of the collective variable ^(x) is the distance 
between the center x c of the nano-particle (a point kept fixed during the simu- 
lations) and the closest oxygen atom, i.e. &(x) = min/{|x c — x?|}, where xf is 
the coordinate of the z-th oxygen atom. The force acting on the atoms associated 
to this collective variable cannot be straightforwardly evaluated since &(x) is a 
non-analytical function of x and, therefore, there is no way toproceed through 
the direct calculation of the V£/2(^(x) -^*) 2 term of Eq. [I[ We replace the 
above definition of the collective variable by a smooth explicit approximation of 
it that, in a proper limit, converges to its exact definition. This smooth explicit 
approximation is obtained in two steps: (i) first we obtain an explicit expression 
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of min/{|x c — xP|}and (ii) then we introduce a smooth approximation to this ex- 
pression. The first step consists in recognizing the following identity: 



min{|x c -x?|} = £ 



*?\) 



(5) 



where &{x) is the Heaviside step function 



0(x) 



x<0 
x>0 



(6) 



x? 



Let / be the O atom closest to the center of the nano-particle, then YljL t 0(\x c — 

( — |x c — x?|) = 8u , where 8u is the Kronecker symbol. This implies that the 

result of the r.h.s. of Eq. [5 is |x c —xf\, i.e. the distance from x c of the closest O 
atom. Eq. [5] is, therefore, trie explicit expression of the collective variable 
A smooth approximation to ^(x) can be obtained by replacing the Heaviside step 
function by a sigmoid function. We use a sigmoid function expressed in terms of 
the Fermi function: 



5(0 = 1- 



1 +exp[Ar] 



(7) 



where A is the parameter controlling its smoothness. In our simulations A has been 
chosen such that the sigmoid function goes from 0.95 to 0.05 in one atomic layer 
(rsj 2 A). A consequence of this replacement is that the size of the nano-particle is 
now defined as a weighted average of the distance of one atomic layer of oxygen 
atoms from the centre of the nano-particle. 



3.2 The ^( x ) collective variable to monitor the phase of the nanoparticle 

We compute the free energy of a Si nano-particle embedded in silica as a func- 
tion of its degree of order, as measured by the bond-orientational order parameter 
(i?6( x )) introduced by Steinhardt et al. ifTTl for bulk systems. In this work we us a 
modified version of the original definition adapted to the case of confined systems, 
as described below in detail. 

In general, JS^ (x) is defined as 



^6(x)=(^-r £ |^ 6m (x)| 2 ) (8) 
\ m=— 6 ' 

where ^ m (x) is the normalized and weighted sum of q l 6m (x) (defined below) 
which, in bulk systems, reads 

*.(.)- m 
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where N is the number of atoms in the system, Ni is the number of nearest neigh- 
bors of the atom i and m = —6, . . . , 6. In the case of confined systems we limit the 
sum over i to just the atoms belonging to the nanoparticle. Consistently with our 
definition of the size of the nanoparticle, we identify these atoms as those at a dis- 
tance lower than ^* from the center of the nanoparticle (^* is the restrain value 
for the collective variable &(x) used in the RMD-PT simulation). According to 
this definition, the ^m(x) of the nanoparticle is: 

„ , x E^i^L(x)G)(^*-|xf -x c |) 

^6m(x) = U ~ X ^ 6mV ' 1 ^ (10) 

As for the case of the collective variable ^(x), the Heaviside step function is 
replaced by a sigmoid function S(t) = 1 — 1/(1 + exp[Af]). 



The q l 6m {x) function appearing in Eq. 10 is defined according to the following 
expression: 

i ( \ ^J= l Y 6m(xij) 
<?6m( x ) = N (11) 

where Y^ m (xij) is the spherical harmonic function of degree 6 and component m 
computed on the solid angle Xjj formed by the distance vector Xij = x; — x; and 
the reference system. The sum runs over the Ni nearest neighbors of the atom 
i. The sum over the m component in Eq. [8] makes the collective variable ^(x) 
rotationally invariant, i.e. independent on the orientation of the reference system. 

When the system is an ideal crystal and the temperature is K, the environ- 
ment of all the atoms is the same and, therefore, J?6 ( x ) is maximum as there is not 
interference among the q l 6m (x). The value of the ^(x) depends on the structure 
of the crystal. For bulk Si JS^x) =0.63. On the contrary, in a perfectly disor- 
dered system the orientation of bonds is random and, therefore, there is complete 
interference among the q l 6m (x), and ^(x) is zero - However, even when the sys- 
tem is at finite temperature and its size is finite, this order parameter is still able 
to distinguish between a disordered and a crystalline phase, being its fluctuation 
(°^ 6 ~ 0.005) much smaller than the difference A JS^ between the two states (typ- 
ically A£? 6 >0.1, see Sec.|4]). 



3.3 Simulation setup 

In the present investigation, the restrained MD is governed by the superposition of 
the Billeter et al. environment-dependent force field 1 18 ] and the restraining po- 
tential k/2(£! 6 (x) - 3D 2 +k' /2(&(x)-&*) 2 . k and k' are the parameters con- 
trolling the degree of biasing and must be large enough such that along the MD the 
values of JS^x) and &(x) oscillate about the target values £l\ and respec- 
tively. As for the Billeter et al. potential, its reliability in modeling equilibrium 
and dynamical properties of Si nanoparticles embedded in silica, of the Si/a-Si02 
interface and Si nanowires has been already well established B18lll9ll20ll21ll22ll . 

The samples are prepared by thermally annealing a periodically-repeated amor- 
phous silica system and embedding Si nanograins (extracted from a well equili- 
brates either amorphous or crystalline bulk). The amorphous silica sample is pre- 
pared through the quenching-from-the-melt procedure, that is by cooling down 
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very slowly a high temperature Si02 melt. The total system contains from ~ 6000 
to ~ 12000 particles, corresponding to a nano-particle radius varying in the range 
1—2 nm. The samples are first thermalized at 300 K in order to release possible 
stress at the Si/silica interface. Typically, during such a thermalization step, the 
nanoparticles slightly shrink. After this initial step, we impose the restraint on the 
size of the nanoparticles and thermalize the samples at the various target temper- 
atures. After this treatment the samples are ready for the restrained simulations 
described above. 

In order to verify possible artifacts due to finite- size effects, we repeated the 
calculation of the mean force at few selected values of £H\ and ^* on samples at an 
increasing size of the silica matrix. We did not observe any significant difference 
in them (the differences were within the statistical error). This demonstrates that 
there are no finite- size effects in our free energy calculations. 

As for the Parallel Tempering temperatures, we used the following seven val- 
ues: 300 K, 500 K, 750 K, 1000 K, 1250 K, 1500 K, 1750 K. 



4 Results and discussion 

In order to compute the phase diagram of the embedded Si nanoparticle, the first 
point is to find minima in the free energy F(^\8^), where F(J2g\&) denotes the 
free energy as a function of ^6 for a Si nanoparticle of size 8%. In the present 
case, since we deal only with two collective variables, and since we know from 
experimental and theoretical facts what is the range of the values of the collective 
variables (see below), we can directly proceed to the reconstruction of the free 
energy in a given interval of both and 8%. We compute the mean force at 
a set of points in the range JS^ G [0,0.65] and 8? E [0.8, 1.8] nm. The rationale 
for the upper limit of the i?6 range is that the value of the bond-orientational 
order parameter for an ideal Si crystal is = 0.63 and, since from experiments 
and previous calculations it is known that Si crystalline nanoparticles assume a 
structure with a (distorted) diamond-like core and a disordered periphery 1123 11241 
[25ll26lh we expect the JS^ of crystalline nanoparticles be lower than this limit. The 
samples created according to the protocol described in Sec. |3.3| confirm that the 
i?6 of crystalline nanoparticles is lower than this upper limit. However, after a 
preliminary scanning of V^ 6 F(^\8?) that allowed to identify the region of JS^ 
containing the minima in the above domain, we restricted the calculations to a 
smaller range: [0.04,0.285], [0.07,0.35] and [0.03,0.38] for the nanoparticles of 
radius 0.8, 1.3 and 1.8 nm, respectively. As for the size of nanoparticles, the same 
experiments mentioned above report unusual phase transitions (disorder-to-order 
with growing T) for nanoparticles in the range [0.5, 2.0] nm. We decided therefore 
to study nanoparticles of three size in this range: 0.8, 1.3 and 1.8 nm. 

We initially performed our calculations without using the combined RMD-PT 



approach (i.e. purely restrained MD as described in Sec. 2.1 ). An example of the 
F{&(\8%) function resulting from these calculation is shown in Fig. [3] where the 
details of the simulation protocol are as follow. We start our simulations from 



one configuration either in the crystalline or amorphous domain (see Sec. |3.3| ), 
set the value of J?6 consistently with this configuration and, after a period of re- 
laxation, compute the mean force at this value of The initial configuration 
for the next value of £?6 is selected randomly from the restrained MD run at the 
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Fig. 3 Free energy vs curves for a nanoparticle of radius ^ = 1 .3 nm at T — 750 K. The free 
energy is calculated along the amorphous to crystalline (O) and the crystalline to amorphous (□) 
paths, as described in the text. The free energy curve for the same nanoparticle size and system 
temperature as obtained from a RMD-PT simulation (A) is also reported for comparison. An 
offset has been applied to this latter curves to improve readability. 



previous value of Of course, a preliminary restrained MD is ran before com- 
puting the mean force, so as to let the system to relax to the new value of i?6. In 
the following, we will refer to the sequence of initial configurations as a "path". 
So, for example, the sequence selected starting from a low J?6 state an d going 
to a high ^6 state will be referred to as Amorphous — >> Crystalline path. The 
opposite path will be called Crystalline — >> Amorphous. Back to the analysis of 
the purely restrained MD results shown in Fig. [3| we notice that the mean force 
computed along the Amorphous —> Crystalline path is different from that com- 
puted going in the backward direction. Moreover, the system presents only one 
minimum: in the amorphous domain when the mean force is computed along the 
Amorphous —> Crystalline path, in the crystalline domain in the other case. To 
tackle this situation, which clearly points out to an hysteresis of the second type 
illustrated in Sec. |2.1[ we repeated the same simulation described above with the 
RMD-PT method and found that (i) the F(£l 6 \M) is independent of the path di- 
rection and (ii) it presents two minima, one at low and one at high £?6 (see Fig. [3}. 

We now move to the analysis of the results of our RMD-PT simulations. In 
Fig. [4] the free energy curves F(£}g\&) of Si nano-particles of size 8% = 0.8 nm, 
^=1.3 nm, and ^=1.8 nm at various temperatures in the overall range 227° C - 
1477°C are shown (the results are presented in Celsius for homogeneity with avail- 
able experimental data). As a first remark, we notice that two metastable states are 
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Amorphous O Cnsmlliiic Amorphoi Q 6 Crystalline Amorphous Q 6 Crystalline 



Fig. 4 Free energy vs £?6 curves for nano-particles with radius 0.8 nm, 1.3 nm and 1.8 nm. The 
curves are shifted to improve readability. 



present at all temperatures and sizes. In general, to high values of J?6 correspond 
crystalline states while to low values of the same order parameter correspond dis- 
ordered (amorphous) states. However, especially for the smallest nanoparticle, the 
difference between the value of J2e corresponding to the two metastable states 
is small. Therefore, the identification of the phase corresponding to a given state 
requires a further investigation of the corresponding structure. We performed this 
analysis by computing the Si-Si partial pair correlation function g(R) for the atoms 
belonging to the nanoparticle (|xf — x c | < ^*) averaging over configurations cor- 
responding to the two metastable states. In the left-most panel of Fig. [5] we report 
the g(R) of the low (bottom panel) and high (top panel) JS^ metastable states for 
the 8% = 0.8 nm nanoparticle at various temperatures. For the sake of comparison, 
we also show the g(R) of bulk amorphous and crystalline states at T = 627°C. 
For the high £?6 metastable state, we notice that even at the highest T the g(R) is 
characterized by three peaks in the range < R < 5 A, one of which is sometimes 
barely visible. These peaks correspond to the bulk-like first, second, and third 
neighboring shell, respectively. They broaden and reduce in intensity by increas- 
ing the temperature; nevertheless, they are still well visible even at T = 1227°C. 
In general, even at low T, these peaks are broader than the corresponding bulk 
crystalline ones. As for the low 82§ metastable state, we notice that in the same R 
range analyzed for the high £?6 case there are only two peaks. The first one, sharp 
and intense, corresponds to the nearest neighbor Si-Si pairs. The second one, very 
broad and small, is usually assigned to the second neighbor pairs, which in amor- 
phous system are distributed over a broad R range. There is no other peak in the 
< R < 5 A domain. Comparing the g(R) of the low £?6 metastable state of this 
particle with amorphous bulk Si we notice that there is a one to one correspon- 
dence between the equivalent peaks in the two system. Similar results are found 
for the 8% = 1 .3 nm and ^=1.8 nm nanoparticles (see the central and right-most 
panels of Figs. [5] respectively). On the basis of this analysis, we identified the high 
metastable state of all the nanoparticles at all temperatures as crystalline and 
the one at low JS^ as amorphous. 

The above conclusions are confirmed by the visual inspection of the structure 
of the nanoparticles in the low and high domains. In Fig. [6] we show few 
snapshots collected along our RMD-PT simulations at JS^ values corresponding 
to the minima of the free energy at low and high temperatures. It appears evident 
that the nanoparticles in the high J?6 domain have a crystal-like core in which the 
tetrahedral orientation of the Si-Si bonds is preserved. This core is surrounded by a 
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Fig. 5 Si-Si pair correlation function (g(r)) of the low and high metastable states at various 
T. Bulk crystalline and amorphous g(r) are also reported for comparison. 



disordered shell. This finding is consistent with previous results (26). Instead, the 
structure of the nanoparticles corresponding to the low JS^ free energy minimum 
is completely disordered. 

Having identified the nature of the metastable states we turn to the study of 
the phase diagram of Si nanoparticles with size and temperature of the nanopar- 
ticle. Our simulations (see Fig. [4]) provide the following qualitative sharp picture: 
for small nano-particles {£% = 0.8 — 1.3 nm) at low temperature (T < 727°C) the 
most stable phase is the amorphous one. The crystalline phase becomes equiprob- 
able (3% = 0.8 nm) or more probable (3% = 1.3 nm) at higher temperatures. On 
the contrary, for larger particles (3% > 1.8 nm) the behavior is bulk-like: at low 
temperatures (T < 977° C) the crystalline phase is the most stable one while at 
higher temperatures the amorphous phase is preferred . 

Let us quickly summarize some experimental findings that we will use for 
comparison. Results from Energy Filtered Transmission Electron Microscopy (EFTEM) 
and Dark-Field Transmission Electron Microscopy (DFTEM) on Si-rich SiO* 
show that Si nanoparticles start to form at 1000°C 128U29ll24U25l[30l[27ll . At this 
temperature they are all amorphous, while at 1 100°C about one third become crys- 
talline. By further increasing the temperature by 50°C, the fraction of crystalline 
nanoparticles rises up to 60%, while the average size of the nano-particles and 
the distribution of their size remains almost unchanged. Finally, at 1250°C 100% 
of nano-particles are crystalline. At this temperature the average size is slightly 
increased, but the particle size distribution still overlaps with the distributions at 
1 100°C and 1 150°C. Our results provide the following interpretation of the these 
experimental results. At low temperatures the nanoparticles are small (as shown by 
experiments) and amorphous due to the found inversion of stability for nanopar- 
ticle of this size. At moderately higher temperatures, when size distribution of 
the nanoparticles is essentially unchanged, the larger nanoparticles in the sample 
transform from amorphous to crystalline, as the amorphous-to-crystalline transi- 
tion temperature is lower for these particles. By further increasing the temperature 
the average size of the nano-particles increases significantly and almost all the 
nanoparticles transform into crystalline as at this size and temperature this is the 
most stable phase. The remaining small nanoparticles are crystalline too due to 
the inversion of stability. In other words, our results bring to the conclusion that 
the observed crystallization of the nanoparticles with the increase of the temper- 
ature is due to the interplay of the effect of the temperature on their size and the 
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TZ* = 0.8 nm 



T = 27°C T = 977°C 



1.8 nm 




T=27°C T = 977°C 



T = 227°C T = 977°C 



Fig. 6 (color online) Snapshots of the nanoparticle of various size in the amorphous and crys- 
talline metastable states at low and high temperature. Light (yellow) and dark (red) spheres 
represent the Si and O atoms, respectivelly. 



inversion of the relative stability of the amorphous and crystalline phase with the 
temperature for small nanoparticles. 



5 Conslusions 

We have introduced a combined restrained MD/Parallel Tempering approach to 
compute the difference in free energy as a function of a (a set of) collective vari- 
able between two metastable states in presence of unknown slow degrees of free- 
dom. We have applied the RMD-PT approach to the study of the phase diagram of 
Si nanoparticles embedded in amorphous silica as a function of their size and the 
system temperature. We found that, at variance with bulk Si, at low T the small 
nanoparticles are amorphous, undergoing an amorphous-to-crystalline phase tran- 
sition at higher T. Large nanoparticles instead recover, as expected, the bulk- like 
behavior; crystalline at low T and amorphous at higher T. We also investigated 
the structure of the nanoparticles, finding that, in agreement with previous works 
ll26lL the crystalline nanoparticles are formed by an ordered core surrounded by a 
disordered periphery. 
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